##Make tables and graphs that haven't been made in other programs

rm(list=ls())

source('param_general.R')
age06Max <- ageCut06

source("./SwitchCosts/stata_results_switch.R")

resultsDatFile <- honoreResList


##Title of different specifications
mainResultNm <- 'Main'
honoreRobustNm <- c('Normal_LD','Age_18','Same_Doc','Commercial','MedicaidFFS','NormalVaginal_LD','NormalCS_LD','Movers','ComplexFirst')
honoreRobustNm3 <- c('age18','normalLD','commercial','medicaidFFS','movers','sameClinician')
naiveRegNm <- c('Naive:All','honore')
          
##Load Data
load(paste0(datDir, 'estimDat_minPpl_all_',minZipHospPpl,'.RData'))
useDatFinal <- readRDS(paste0(datDir,'ld_women_process_all.rds', sep = ''))
load(paste0(datDir,'chambDat_minPpl_',"all",'_',minZipHospPpl,'.RData', sep = '')) 
chambDat <- chambDat[age06 <= ageCut06 & noMove3==1 & hosp_id_alt_1_0 != hosp_id_alt_1_1 ,]

##Descriptive Stats (text of data section)
estimDat[,chosenPrevInChoiceSet := sum(chosenPrev)>0, by = 'admission']
message('Number of admission for a women with one child, where the same hospital was chosen as for the previous labor where previous hospital is in choice set',estimDat[birthNum>1 & age06<=ageCut06 &obs_choice == 1 & hosp_id_alt !=0 & chosenPrevInChoiceSet ==TRUE,percent(mean(chosenPrev))])
message('Same, but for normal LD: ',estimDat[birthNum>1 & age06<=ageCut06 &obs_choice == 1 & hosp_id_alt !=0 & normalLD_indic == TRUE & chosenPrevInChoiceSet,percent(mean(chosenPrev))])
message('Same, but for complicated LD: ',estimDat[birthNum>1 & age06<=ageCut06 &obs_choice == 1 & hosp_id_alt !=0 & normalLD_indic != TRUE & chosenPrevInChoiceSet,percent(mean(chosenPrev))])
message('Number of admission for a women with one child, where the same hospital was chosen as for the previous labor (irrespective of whether the previous hospital was in the choice set)',estimDat[birthNum>1 & age06<=ageCut06 &obs_choice == 1 & hosp_id_alt !=0 ,percent(mean(chosenPrev))])

prevHospByBirth <- estimDat[birthNum>1 & age06<=ageCut06 &obs_choice == 1 & hosp_id_alt !=0 & chosenPrevInChoiceSet == 1 ,
                            list(meanShr = mean(chosenPrev), sdDev = sd(chosenPrev), .N),by=birthNum][order(birthNum)]
prevHospByBirth[,stdErr := sdDev/sqrt(N)]
prevHospByBirth[,lbCI := meanShr - stdErr*1.96]
prevHospByBirth[,ubCI := meanShr + stdErr*1.96]
print(prevHospByBirth)


message('Number of women meeting age cutoff with three births ',nrow(estimDat[birthNum == 3 &age06<=ageCut06 & obs_choice ==1]))
message('Share of women meeting age cutoff with three births ',percent(nrow(estimDat[birthNum == 3 &age06<=ageCut06 & obs_choice ==1])/length(unique(estimDat[age06<=ageCut06, maskssn]))))

message('Of those women who switched hospitals between first two births, and both were available both times ',nrow(estimDat[birthNum == 3 &age06<=ageCut06 & obs_choice ==1 & hospSwitchFirst2 == 1]))
message('Of those, women who did not move bw two and three ',nrow(estimDat[birthNum == 3 &age06<=ageCut06 & obs_choice ==1 & noMove3==1 & hospSwitchFirst2 == 1]))

message('Of those, women who are ABA ',nrow(estimDat[birthNum == 3 &age06<=ageCut06 & obs_choice ==1 & firstHosp3==1 & hospSwitchFirst2==1 & noMove3 == 1]))
message('Number of those women who are ABB ',nrow(estimDat[birthNum == 3 &age06<=ageCut06 & obs_choice ==1 &hospSwitchFirst2==1& secondHosp3==1 & noMove3==1]))
message('This is the set used to estimate switching costs')

message('Of that set, this number did not move at all ',nrow(estimDat[birthNum == 3 &age06<=ageCut06 & obs_choice ==1 & noMove3==1 & noMove2==1 & hospSwitchFirst2 == 1]))
message('Of those, women who are ABA ',nrow(estimDat[birthNum == 3 &age06<=ageCut06 & obs_choice ==1 & firstHosp3==1  & noMove3 == 1 & noMove2==1 & hospSwitchFirst2 == 1]))
message('Of those, women who are ABB ',nrow(estimDat[birthNum == 3 &age06<=ageCut06 & obs_choice ==1 & secondHosp3==1  & noMove3 == 1 & noMove2==1 & hospSwitchFirst2 == 1]))


##Summary Stats (tab:Summary Statistics)
cbsaXWalkDat <- read.dta(cbsaDatPath)
allPat <- unique(estimDat[,maskssn])
ageRestPat <- unique(estimDat[age06 <= age06Max,maskssn])
threeBirthPat <- unique(estimDat[age06 <= age06Max & birthNum == 3,maskssn])
chambPat <- unique(chambDat[,maskssn])

useDatFinalMerge <- data.table(merge(data.frame(useDatFinal),cbsaXWalkDat[,c('zcta5','memi')], by.x = 'pat_zip', by.y = 'zcta5',all.x = TRUE))
setkey(useDatFinalMerge, maskssn)

dataList <- list(All = useDatFinalMerge[J(allPat)],
                 AgeRest = useDatFinalMerge[J(ageRestPat)],
                 Honore = useDatFinalMerge[J(chambPat)])

summaryStat <- ldply(dataList, function(x) x[,list(Age =round(mean(age),2),
                                                   White = round(mean(white),2),
                                                   Black = round(mean(black),2),
                                                   Hispanic = round(mean(hispanic),2),
                                                   Medicaid = round(mean(payer %like% 'Medicaid'),2),
                                                   Commercial = round(mean(payer %like% 'Commercial'),2),
                                                   Metropolitan = round(mean(memi == 'Metropolitan' & !is.na(memi)),2),
                                                   `Normal Birth` = round(mean(normalLD_indic),2),
                                                   `Number of Births` = round(length(maskssn)),
                                                   `Number of Women` = round(length(unique(maskssn))))])

summaryStatTypNames <- summaryStat$.id
summaryStat$.id <- NULL
summaryStatVarNames <- colnames(summaryStat)

summaryStatFmt <- t(as.data.table(c(lapply(summaryStat[,1:8],formatC,digits = 2, format = 'f'),
                                  lapply(summaryStat[,9:10],formatC,digits = 0, format = 'f',big.mark = ','))))
colnames(summaryStatFmt) <- summaryStatTypNames

summaryTable <- xtable(summaryStatFmt)
print(summaryTable,floating=FALSE)
align(summaryTable) <- c('r',rep('c',ncol(summaryTable)))
print(summaryTable, file = paste0(resultsDir,'Switch/summaryStatistics.tex'),floating = FALSE)


##Compare women (who meet age restriction) who switched hospital after first birth to women who didn't
message('Table with share of women with complications after first birth, by hospital switching behavior')
print(estimDat[age06<=ageCut06 & obs_choice == 1 & birthNum == 2 & hosp_id_alt!=0,  list(.N,`Normal First Birth` = mean(normalLD_hist),se = sd(normalLD_hist)/sqrt(.N)),by = list(`Different Hospital For Second Birth` = chosenPrev)])


##Estimation results (tables and graphs)
honoreFeResults<- ldply(honoreResList$distAndSwitch, function (x) c(coef(summary(x))[,1:2],N = nobs(x)))
names(honoreFeResults)[names(honoreFeResults) == '.id'] <- "Type"

setnames(honoreFeResults,c("V1","V2","V3","V4"),c("Estimate","Estimate Distance","Std. Error", "Std. Error Distance"))


honoreFeResultsRatio  <- ldply(honoreRatioList$distAndSwitch, function (x) x)
names(honoreFeResultsRatio)[names(honoreFeResultsRatio) == '.id'] <- "Type"
setnames(honoreFeResultsRatio,c("SE"),c("Std. Error"))
honoreFeResultsRatio = as.data.table(cbind(honoreFeResultsRatio,N = honoreFeResults$N))
honoreFeResultsRatio[,Estimate := -1 * Estimate]
honoreFeResultsRatio[,V5 := 0]
honoreFeResultsRatio[,V6 := 0]
setnames(honoreFeResultsRatio,c("V5","V6"),c("Estimate Distance", "Std. Error Distance"))

noColumns <- length(colnames(stataRes))
resultsNaive <- as.data.table(t(rbind(colnames(stataRes)[2:noColumns],
                                      t(as.numeric(stataRes[5,2:noColumns])),t(as.numeric(stataRes[3,2:noColumns])),
                                      t(as.numeric(stataRes[6,2:noColumns])),t(as.numeric(stataRes[4,2:noColumns])),
                                      t(as.numeric(stataRes[7,2:noColumns])))))

noColumnsFeas <- length(colnames(stataFeasibleRes))
resultsFeasibleNaive <- as.data.table(t(rbind(colnames(stataFeasibleRes)[2:noColumnsFeas],
                                      t(as.numeric(stataFeasibleRes[3,2:noColumnsFeas])),NA,
                                              t(as.numeric(stataFeasibleRes[4,2:noColumnsFeas])),NA,
                                      t(as.numeric(stataFeasibleRes[5,2:noColumnsFeas])))))
resultsFeasibleNaive[V1 == 'honore',V1:='honoreFeas']
resultsFeasibleNaive[V1 == 'nonFirstBirth',V1:='nonFirstBirthFeas']

resultsNaive <- rbind(resultsNaive,resultsFeasibleNaive)

setnames(resultsNaive,c("V1","V2","V3",'V4',"V5","V6"),c("Type","Estimate","Estimate Distance","Std. Error", "Std. Error Distance",'N'))
resultsNaive[Type == "base", Type := 'Naive:All']
resultsNaive = rbind(resultsNaive,subset(honoreFeResults, Type == 'Main'))


resultsNaiveRatio = as.data.table(t(rbind(colnames(stataRes)[2:noColumns],
                                     t(as.numeric(stataRes[8,2:noColumns])),t(as.numeric(stataRes[9,2:noColumns])),
                                     t(as.numeric(stataRes[6,2:noColumns])))))

setnames(resultsNaiveRatio,c("V1","V2","V3",'V4'),c("Type","Estimate","Std. Error",'N'))
resultsNaiveRatio[Type == "base", Type := 'Naive:All']
resultsNaiveRatio[,Estimate := -1 * as.numeric(Estimate)]
resultsNaiveRatio[,V5 := 0]
resultsNaiveRatio[,V6 := 0]
setnames(resultsNaiveRatio,c("V5","V6"),c("Estimate Distance", "Std. Error Distance"))
resultsNaiveRatio = rbind(resultsNaiveRatio,subset(honoreFeResultsRatio, Type == 'Main'))

resultsList <- list(honore = data.table(rbind(honoreFeResults, subset(honoreFeResults, Type == 'Naive:All'))),
                    honoreRatio = data.table(rbind(honoreFeResultsRatio, subset(honoreFeResultsRatio, Type == 'Naive:All'))),
                       naive = as.data.table(resultsNaive),
                       naiveRatio = as.data.table(resultsNaiveRatio))



lapply(names(resultsList),function(x) resultsList[[x]][,regClass := x])
resultsListComb <- rbindlist(resultsList,use.names = TRUE)
for (col in c('Estimate','Std. Error','N',"Estimate Distance","Std. Error Distance")) resultsListComb[,(col) := as.numeric(get(col))]
resultsListComb[,':='(lbCI = Estimate - 1.96 * `Std. Error`, ubCI = Estimate + 1.96 * `Std. Error`)]

resultsListComb[Type == 'Main', titleCol := "Fixed Effect"]
resultsListComb[Type == 'Naive:All', titleCol := "Lagged Dep Var \n (Full)"]
resultsListComb[Type == 'honore', titleCol := "Lagged Dep Var \n (Restricted)"]
resultsListComb[Type == 'Age_18', titleCol := "Age Cutoff \n at 18"]
resultsListComb[Type == 'Normal_LD'|Type == 'normalLD', titleCol := "Normal Labor \n and Delivery"]
resultsListComb[Type == 'Same_Doc' | Type == 'sameClinician', titleCol := "Same Clinician"]
resultsListComb[Type == 'Commercial', titleCol :="Commercial Insurance"]
resultsListComb[Type == 'Medicaid', titleCol :="Medicaid Insurance"]
resultsListComb[Type == 'MedicaidFFS'| Type == 'medicaidFFS', titleCol :="Medicaid FFS"]
resultsListComb[Type == 'ComplexFirst', titleCol :="Complicated 1st Birth"]
resultsListComb[Type == 'Movers'|Type == 'movers' , titleCol :="Moved Residences"]
resultsListComb[Type == 'jacksonvilleMSA', titleCol :="Jacksonville MSA"]
resultsListComb[Type == 'NormalCS_LD' , titleCol :="Normal C Section"]
resultsListComb[Type == 'NormalVaginal_LD' , titleCol :="Normal Vaginal"]
resultsListComb[Type == 'age18', titleCol :="Age Cutoff \n at 18"]
resultsListComb[Type == 'commercial', titleCol :="Commercial Insurance"]
resultsListComb[Type == 'medicaid', titleCol :="Medicaid Insurance"]
resultsListComb[Type == 'jacksonville', titleCol :="Jacksonville MSA"]
resultsListComb[Type == 'medicaid', titleCol :="Medicaid"]
resultsListComb[Type == 'normalCS', titleCol :="Normal C Section"]
resultsListComb[Type == 'normalVag', titleCol :="Normal Vaginal"]
resultsListComb[Type == 'secondBirth', titleCol :="Second Births Only"]
resultsListComb[Type == 'honoreFeas', titleCol :="Group Level Heterogeneity \n (Restricted)"]



resultsListComb[,isMain := Type == mainResultNm]
resultsListComb[,titleColTbl := str_replace(titleCol,' \n ',' ')]
          
                               
makeGraphFunc <- function(dat,isRobust){
    
    filePref <- dat[,list(regClass)][1]
    reverseDum = -1
  #reverseDum <- ifelse(names(datList)=='chamb', -1, 1)
  
  if(isRobust == 1) {
    specificationSet = c(naiveRegNm,mainResultNm)
    fileName = paste0(filePref,"SC")
  } 
  if(isRobust == 2) {
    specificationSet = honoreRobustNm
    fileName = paste0(filePref,"SC","Robust")
  }
  if(isRobust == 3) {
    specificationSet = honoreRobustNm3
    fileName = paste0(filePref,"SC","Robust2")
  }
    pltDat <- dat[Type %in% specificationSet,]
  graph <- ggplot(pltDat, aes(x = titleCol, y = Estimate, ymin =  lbCI,
                                                            ymax = ubCI, color = isMain
                              ), environment = environment()) +
    geom_pointrange(size = .75) + scale_x_discrete() +
    guides(color = FALSE)+ scale_color_manual(values= c("blue","red")) + coord_flip()  + theme_bw() +
    labs(y = "Estimate", x = "" ) +
    theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 14)) +
    theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 14)) +
    theme_bw() + theme(panel.grid.major = element_blank(),
                       panel.grid.minor = element_blank() ) + scale_y_continuous(breaks = pretty_breaks(n= 6))
  
  if(isRobust > 1) {
    graph = graph +     geom_hline(yintercept = unlist(dat[dat$Type == mainResultNm,'Estimate',with=FALSE]), color = 'blue') +
      geom_hline(yintercept = unlist(dat[dat$Type == mainResultNm,c('lbCI','ubCI'),with=FALSE]), color = 'blue',linetype = 'longdash') +
      geom_hline(yintercept = unlist(dat[dat$Type == 'Naive:All','Estimate',with=FALSE]), color = 'red') +
      geom_hline(yintercept = unlist( dat[dat$Type == 'Naive:All',c('lbCI','ubCI'),with=FALSE]), color = 'red',linetype = 'longdash') + 
      scale_color_manual(values= c("black"))
  }
  
  ggsave(paste0(resultsDir,"/Switch/",fileName,'Paper.png'),graph,width = widthPaper,
         height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")
  ggsave(paste0(resultsDir,"/Switch/",fileName,'Pres.png'),graph,width = widthPres,
         height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo") 
}


makeGraphFuncBW <- function(dat,isRobust){
  
  filePref <- dat[,list(regClass)][1]
  reverseDum = -1
  #reverseDum <- ifelse(names(datList)=='chamb', -1, 1)
  
  if(isRobust == 1) {
    specificationSet = c(naiveRegNm,mainResultNm)
    fileName = paste0(filePref,"SC")
  } 
  if(isRobust == 2) {
    specificationSet = honoreRobustNm
    fileName = paste0(filePref,"SC","Robust")
  }
  if(isRobust == 3) {
    specificationSet = honoreRobustNm3
    fileName = paste0(filePref,"SC","Robust2")
  }
  pltDat <- dat[Type %in% specificationSet,]
  graph <- ggplot(pltDat, aes(x = titleCol, y = Estimate, ymin =  lbCI,
                              ymax = ubCI, color = isMain
  ), environment = environment()) +
    geom_pointrange(size = .75) + scale_x_discrete() +
    guides(color = FALSE)+ scale_color_manual(values= c("black","darkgrey")) + coord_flip()  + theme_bw() +
    labs(y = "Estimate", x = "" ) +
    theme(axis.title.x = element_text(face = "bold", size = 14, vjust = .1), axis.text.x = element_text(size = 14)) +
    theme(axis.title.y = element_text(face = "bold", size = 14, vjust = .1), axis.text.y = element_text(size = 14)) +
    theme_bw() + theme(panel.grid.major = element_blank(),
                       panel.grid.minor = element_blank() ) + scale_y_continuous(breaks = pretty_breaks(n= 6))
  
  if(isRobust > 1) {
    graph = graph +     geom_hline(yintercept = unlist(dat[dat$Type == mainResultNm,'Estimate',with=FALSE]), color = 'black') +
      geom_hline(yintercept = unlist(dat[dat$Type == mainResultNm,c('lbCI','ubCI'),with=FALSE]), color = 'black',linetype = 'longdash') +
      geom_hline(yintercept = unlist(dat[dat$Type == 'Naive:All','Estimate',with=FALSE]), color = 'darkgrey') +
      geom_hline(yintercept = unlist( dat[dat$Type == 'Naive:All',c('lbCI','ubCI'),with=FALSE]), color = 'darkgrey',linetype = 'longdash') + 
      scale_color_manual(values= c("black"))
  }
  
  ggsave(paste0(resultsDir,"/Switch/",fileName,'PaperBW.png'),graph,width = widthPaper,
         height = heightPaper, dpi = dpiPlot, units = unitsPaper, type = "cairo")
  ggsave(paste0(resultsDir,"/Switch/",fileName,'PresBW.png'),graph,width = widthPres,
         height = heightPres, dpi = dpiPlot, units = unitsPres, type = "cairo") 
}



#resultsLogList['chamb'] c(naiveRegNm,mainResultNm)
#specificationNamesList[[1]]

makeGraphFunc(resultsListComb[regClass == 'honore' | regClass == 'naive'],isRobust = 1)
makeGraphFunc(resultsListComb[regClass == 'honoreRatio' | regClass == 'naiveRatio'],isRobust = 1)
makeGraphFunc(resultsListComb[regClass == 'honore' | regClass == 'naive'],isRobust = 2)
makeGraphFunc(resultsListComb[regClass == 'honoreRatio' | regClass == 'naiveRatio'],isRobust = 2)
makeGraphFunc(resultsListComb[regClass == 'naive'],isRobust = 3)
makeGraphFunc(resultsListComb[regClass == 'naiveRatio'],isRobust = 3)



makeGraphFuncBW(resultsListComb[regClass == 'honore' | regClass == 'naive'],isRobust = 1)
makeGraphFuncBW(resultsListComb[regClass == 'honoreRatio' | regClass == 'naiveRatio'],isRobust = 1)
makeGraphFuncBW(resultsListComb[regClass == 'honore' | regClass == 'naive'],isRobust = 2)
makeGraphFuncBW(resultsListComb[regClass == 'honoreRatio' | regClass == 'naiveRatio'],isRobust = 2)
makeGraphFuncBW(resultsListComb[regClass == 'naive'],isRobust = 3)
makeGraphFuncBW(resultsListComb[regClass == 'naiveRatio'],isRobust = 3)

##Share of women going to previous hospital under different naive estimates
sink(paste0(resultsDir,"Switch/baselineParamEst.tex"))
regResult <- resultsListComb[(regClass == 'naive' & (Type == "baseNoSC" | Type == "Naive:All"))]
cat("\\begin{tabular}{@{\\extracolsep{3pt}}p{3cm}cc}\n")
cat("\\\\[-1.8ex]\\hline \\\\[-1.8ex]\n")
cat(sprintf(" & %s & %s \\\\\n",'Standard Logit','Standard with Lagged Dep Var' ))
cat("\\hline \\hline \\\\[-1.8ex]\n")
cat(sprintf("Travel Time & %3.2f  & %3.2f\\\\\n",regResult[Type=='baseNoSC',`Estimate Distance`],regResult[Type=='Naive:All',`Estimate Distance`]))
cat(sprintf(" & (%3.2f) & (%3.2f)\\\\\\\\\n",regResult[Type=='baseNoSC',`Std. Error Distance`],regResult[Type=='Naive:All',`Std. Error Distance`]))
cat(sprintf("Previous Choice  & & %3.2f \\\\\n",regResult[Type=='Naive:All',Estimate]))
cat(sprintf(" &  & (%3.2f) \\\\\\\\\n",regResult[Type=='Naive:All',`Std. Error`]))
cat(sprintf("Predicted Share to Previous & %3.2f\\%% &  %3.2f\\%% \\\\\n",round(as.numeric(stataRes[stataRes$X=='shrPrevHosp','baseNoSC']),4)*100 ,round(as.numeric(stataRes[stataRes$X=='shrPrevHosp','base']),4)*100))
cat("\\hline\n")
cat(sprintf("N & %s & %s \\\\\n",
            prettyNum(regResult[Type=='baseNoSC',`N`],big.mark = ','),
            prettyNum(regResult[Type=='Naive:All',`N`],big.mark = ',')))
cat("\\hline \\\\[-1.8ex]\n")
cat("\\end{tabular}")
sink()


##Compute pr of going to previous choice under different assumptions
##Zip code level Xis are computed on the basis of first births
load(paste0(datDir,"deltaDat_minPpl_all_5.RData"))
aggDeltDatWide <- dcast.data.table(aggDeltDat,hosp_id_alt~popSet+regTyp,value.var = 'aggXi')
aggDeltDatWide[,hosp_id_alt := as.integer(hosp_id_alt)]
setnames(aggDeltDatWide,setdiff(names(aggDeltDatWide),'hosp_id_alt'),paste0('aggxi_',setdiff(names(aggDeltDatWide),'hosp_id_alt')))
deltDat = grpShrDatXi[is.finite(xiHat),]
setnames(deltDat,"nGrp","nZip")
deltDatWide <- dcast.data.table(deltDat[][,.id := paste0('xi_',.id)],hosp_id_alt + pat_zip ~ .id, value.var = 'xiHat')
setkey(deltDatWide,hosp_id_alt,pat_zip)
setkey(estimDat,hosp_id_alt,pat_zip)
estimXiDat <- merge(estimDat[birthNum > 1& hospSwitchFirst2==1&noMove3==1&birthNum==3,list(time_current,chosenPrev,hosp_id_alt,yearDat,admission,maskssn,pat_zip,obs_choice)],deltDatWide,by = c('hosp_id_alt','pat_zip'),all.x=TRUE)
estimXiDat <- merge(estimXiDat,aggDeltDatWide, by = c('hosp_id_alt'))

estimXiDat[hosp_id_alt != 0,delt_fe_mean := aggxi_honoreSub_fe + time_current * coefDistList$main + coefSwitchList$main * chosenPrev]
estimXiDat[hosp_id_alt != 0,delt_lagDepVar_mean:= aggxi_honoreSub_lagdv + time_current * coefDistList$naiveHonore + coefSwitchList$naiveHonore * chosenPrev]
estimXiDat[hosp_id_alt != 0,delt_noLDV_mean:= aggxi_honoreSub_stlogit+ time_current * coefDistList$noLDVHonore]


estimXiDatMelt <- melt(estimXiDat[,list(hosp_id_alt,pat_zip,chosenPrev,admission,obs_choice,delt_fe_mean,delt_lagDepVar_mean,delt_noLDV_mean)],id.vars = c('hosp_id_alt','pat_zip','chosenPrev','admission','obs_choice'),value.name = 'delta')
estimXiDatMelt[hosp_id_alt == 0, delta := 0]
estimXiDatMelt[,choicePr := exp(delta)/sum(exp(delta)),by = c('admission','variable')]
estimXiDatMelt[chosenPrev == 1 & !is.na(choicePr) ,list(choicePrMean = mean(choicePr,na.rm = TRUE),obsChoiceMean = mean(obs_choice)),by = 'variable']



sink(paste0(resultsDir,"Switch/mainParamEst.tex"))
regResult <- resultsListComb[(regClass == 'honore' | (regClass == 'naive' & Type != "Main")) & Type %in% c(naiveRegNm,mainResultNm)]
regResultRatio <- resultsListComb[(regClass == 'honoreRatio' | (regClass == 'naiveRatio' & Type != "Main")) & Type %in% c(naiveRegNm,mainResultNm)]
cat("\\begin{tabular}{@{\\extracolsep{3pt}}p{3cm}ccc}\n")
cat("\\\\[-1.8ex]\\hline \\\\[-1.8ex]\n")
cat(sprintf(" & %s & %s & %s\\\\\n",regResult[Type=='Naive:All',titleColTbl],regResult[Type=='honore',titleColTbl], regResult[Type=='Main',titleColTbl]))
cat("\\hline \\hline \\\\[-1.8ex]\n")
cat(sprintf("Travel Time & %3.2f  & %3.2f & %3.2f\\\\\n",regResult[Type=='Naive:All',`Estimate Distance`],regResult[Type=='honore',`Estimate Distance`], regResult[Type=='Main',`Estimate Distance`]))
cat(sprintf(" & (%3.2f) & (%3.2f) & (%3.2f)\\\\\\\\\n",regResult[Type=='Naive:All',`Std. Error Distance`],regResult[Type=='honore',`Std. Error Distance`], regResult[Type=='Main',`Std. Error Distance`]))
cat(sprintf("Previous Choice  & %3.2f & %3.2f & %3.2f\\\\\n",regResult[Type=='Naive:All',Estimate],regResult[Type=='honore',Estimate], regResult[Type=='Main',Estimate]))
cat(sprintf(" & (%3.2f) & (%3.2f)  & (%3.2f)\\\\\\\\\n",regResult[Type=='Naive:All',`Std. Error`],regResult[Type=='honore',`Std. Error`], regResult[Type=='Main',`Std. Error`]))
cat(sprintf("Previous Choice  / Travel Time & %3.2f & %3.2f & %3.2f\\\\\n",regResultRatio[Type=='Naive:All',Estimate],regResultRatio[Type=='honore',Estimate], regResultRatio[Type=='Main',Estimate]))
cat(sprintf(" & (%3.2f) & (%3.2f)  & (%3.2f)\\\\\\\\\n",regResultRatio[Type=='Naive:All',`Std. Error`],regResultRatio[Type=='honore',`Std. Error`], regResultRatio[Type=='Main',`Std. Error`]))

cat(sprintf("Hospital/Patient Fixed Effects &  &  & X \\\\\n"))
cat("\\hline\n")
cat(sprintf("N & %s & %s & %s \\\\\n",
            prettyNum(regResult[Type=='Naive:All',`N`],big.mark = ','),
            prettyNum(regResult[Type=='honore',`N`],big.mark = ','),
            prettyNum(regResult[Type=='Main',`N`],big.mark = ',')))
cat("\\hline \\\\[-1.8ex]\n")
cat("\\end{tabular}")
sink()

sink(paste0(resultsDir,"Switch/naiveFeasibleParamEst.tex"))
regResult <- resultsListComb[Type %in% c('honoreFeas','honore')& regClass == 'naive']
cat("\\begin{tabular}{@{\\extracolsep{3pt}}p{3cm}ccc}\n")
cat("\\\\[-1.8ex]\\hline \\\\[-1.8ex]\n")
cat(sprintf(" & %s & %s\\\\\n",regResult[Type=='honore',titleColTbl], regResult[Type=='honoreFeas',titleColTbl]))
cat("\\hline \\hline \\\\[-1.8ex]\n")
cat(sprintf("Previous Choice  &  %3.2f & %3.2f\\\\\n",regResult[Type=='honore',Estimate],regResult[Type=='honoreFeas',Estimate]))
cat(sprintf(" & (%3.2f) & (%3.2f) \\\\\\\\\n",regResult[Type=='honore',`Std. Error`],regResult[Type=='honoreFeas',`Std. Error`]))
cat(sprintf("Previous Choice  / Travel Time & %3.2f & %3.2f \\\\\n",regResultRatio[Type=='honore',Estimate],regResultRatio[Type=='honoreFeas',Estimate]))
cat(sprintf(" & (%3.2f) & (%3.2f)  \\\\\\\\\n",regResultRatio[Type=='honore',`Std. Error`],regResultRatio[Type=='honoreFeas',`Std. Error`], regResultRatio[Type=='Main',`Std. Error`]))
cat("\\hline\n")
cat(sprintf("N & %s & %s  \\\\\n",
            prettyNum(regResult[Type=='honore',`N`],big.mark = ','),
            prettyNum(regResult[Type=='honoreFeas',`N`],big.mark = ',')))
cat("\\hline \\\\[-1.8ex]\n")
cat("\\end{tabular}")
sink()


sink(paste0(resultsDir,"Switch/honoreRobustEst.tex"))
regResult <- resultsListComb[regClass == 'honore'  & Type %in% c(honoreRobustNm)]
regResultRatio <- resultsListComb[regClass == 'honoreRatio'  & Type %in% c(honoreRobustNm)]
cat("\\begin{tabular}{@{\\extracolsep{3pt}}p{4.5cm}cccc}\n")
cat("\\\\[-1.8ex]\\hline \\\\[-1.8ex]\n")
cat(sprintf(" &  Travel Time & Previous Choice & Previous Choice (min) & N \\\\\n"))
cat("\\hline \\hline \\\\[-1.8ex]\n")
for (typ in unique(regResult$Type)){
    cat(sprintf("%s & %3.2f & %3.2f & %3.2f & %s\\\\\n", regResult[Type == typ,titleColTbl],
                regResult[Type == typ,`Estimate Distance`],
                regResult[Type == typ,Estimate],regResultRatio[Type == typ,Estimate],
                prettyNum(regResult[Type == typ,N],big.mark =',')))
    cat(sprintf(" & (%3.2f) & (%3.2f) & (%3.2f) \\\\\\\\\n",regResult[Type == typ, `Std. Error Distance`],
                regResult[Type == typ, `Std. Error`],regResultRatio[Type == typ, `Std. Error`]))
}
cat("\\hline \\\\[-1.8ex]\n")
cat("\\end{tabular}")
sink()

sink(paste0(resultsDir,"Switch/naiveRobustEst.tex"))
regResult <- resultsListComb[regClass == 'naive'  & Type %in% c(honoreRobustNm3)]
regResultRatio <- resultsListComb[regClass == 'naiveRatio'  & Type %in% c(honoreRobustNm3)]
cat("\\begin{tabular}{@{\\extracolsep{3pt}}p{4.5cm}cccc}\n")
cat("\\\\[-1.8ex]\\hline \\\\[-1.8ex]\n")
cat(sprintf(" &  Travel Time & Previous Choice & Previous Choice (min) & N \\\\\n"))
cat("\\hline \\hline \\\\[-1.8ex]\n")
for (typ in unique(regResult$Type)){
  cat(sprintf("%s & %3.2f & %3.2f & %3.2f & %s\\\\\n", regResult[Type == typ,titleColTbl],
              regResult[Type == typ,`Estimate Distance`],
              regResult[Type == typ,Estimate],regResultRatio[Type == typ,Estimate],
              prettyNum(regResult[Type == typ,N],big.mark =',')))
  cat(sprintf(" & (%3.2f) & (%3.2f) & (%3.2f) \\\\\\\\\n",regResult[Type == typ, `Std. Error Distance`],
              regResult[Type == typ, `Std. Error`],regResultRatio[Type == typ, `Std. Error`]))
}
cat("\\hline \\\\[-1.8ex]\n")
cat("\\end{tabular}")
sink()


